Three Dimensional Lattice-Boltzmann Model for Electrodynamics 
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In this paper we develop a 3D Lattice-Boltzmann model that recovers in the continuous limit the 
Maxwell equations for macroscopic mediums. The model can sucessfully reproduces the propagation 
of the electromagnetic waves in dielectric mediums and waveguide, the skin effect, the electrical 
dipole radiation and the electromagnetic response of a resonant cavity. 
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I. INTRODUCTION 

To simulate the electromagnetic fields inside materi- 
als is a fundamental tool in optics and electrodynamics, 
specially when boundary conditions are so complex that 
an analytical solution is impossible. Some examples of 
this phenomena are the design of antennas and the study 
of electrical discharges across media, among many oth- 
ers. For this reason, numerous numerical methods for the 
Maxwell equations have been implemented. One of these 
methods is the FDTD (Finite-difference time-domain) [l| , 
which is based on finite differences and works pretty well, 
but consumes huge ammounts of computational resource. 
Another one is the ray tracing method |2j, that is only 
useful when the electromagnetic signals are of a single 
frequency. In addition, one of the actual problems in 
computational electrodynamics is the simulation of elec- 
tromagnetic pulses in heterogenous media, that is one of 
the main advantages of lattice-Boltzmann methods. 

One of the numerical methods for simulating fluids is 
Lattice Boltzmann (LB) [H, which was developed from 
lattice-gas automata. Lattice Boltzmann simulations are 
performed on regular grids of many cells and a small num- 
ber of velocity vectors per cell, each one associated to a 
set of density distribution functions, which evolve and 
spread together to the ncighbohr cells according to the 
collisional Boltzmann equation. As far as we know, there 
is still not a lattice-Boltzmann model for electrodynam- 
ics, but the effects of Maxwell equations on plasmas has 
been included (in the form of the diffusion equation for 
the magnetic field) as part of LB models for magnetohy- 
drodynamics. The first LB model for studying plasmas 
reproduces the resistive magnetohydrodynamic equations 
and was developed by Chen 0, H| as an extension of the 
Lattice-Gas model developed by Chen and MatthaeusQ 
and Chen, Matthaeus and Klein Q • This LB model uses 
37 velocity vectors per cell on a square lattice and is de- 
veloped for two dimensions. Thereafter, Martinez, Chen 
and Matthaeuss Q decreased the number of velocity vec- 



tors from 37 to 13, which made easier a future 3D exten- 
sion. Some of the first LB models for magnetohydro- 
dynamics in 3D were developed by Bryan R. OsbornQ 
and Breyiannis and Valougeorgi [lfj |. They used 19 vec- 
tors on a cubic lattice for the fluid, plus 7 vectors for 
the magnetic field, which makes a total number of 26 
vectors per cell. By following a different path, Fogaccia, 
Benzi and Romanelli [ll| introduced a 3D LB model for 
simulating turbulent plasmas in the electrostatic limit. 
All these models reproduce the resistive magnetohydro- 
dynamc equations for a single fluid. Nevertheless, in a 
recent work[12], we have introduced a 3D LB model that 
reproduces the two fluids theory for plasmas. This model 
uses a cubic lattice and a set of 39 velocity vectors and 
includes explicitly the Maxwell equations in the vacuum. 

In this paper, we introduce a 3D lattice-Boltzmann 
model that recovers the Maxwell equations in material 
media. The model is based on our previous develop for 
magnetohydrodynamics [12J , but it is able to reproduce 
the behavior of electromagnetic fields inside dielectrics, 
magnets and conductors. The model uses 33 vectors per 
cell and 50 probability density functions, and performs 
well in a wide range of tests. Section fll] describes the 
model, with the evolution rules and the equilibrium ex- 
pressions for the 50 density functions, plus the procedure 
to compute both the charge and current densities and 
the electric and magnetic fields. The Chapman-Enskog 
expansion showing how these rules recover the electrody- 
namic equations is developed in Appendix [Al In order to 
validate the model, we simulate in section ILTIl the reflec- 
tion of electromagnetic pulses on dielectric media, the 
propagation of electromagnetic waves along waveguide, 
the skin effect on conductors, the radiation pattern of an 
oscillating electrical dipole and the natural frequencies of 
a resonant cavity. The main results and conclusions are 
summarized in section ITVl 



II. 3D LATTICE-BOLTZMANN MODEL FOR 
ELECTRODYNAMICS 
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In a simple Lattice-Boltzmann model [3J, the D- 
dimensional space is divided into a regular grid of cells. 
Each cell has Q vectors that links the cell with its 



FIG. 1: Cubic Lattice D3Q13 for modelling the Maxwell equa- 
tions. The arrows represent the velocity vectors vf and the 
electric field vectors e?-, where p indicates the plane of loca- 
tion. 




FIG. 2: Cubic Lattice D3Q7 for simulating the magnetic field. 
The arrows indicate the magnetic vectors b\- . 

neighbors, and each vector is associated to a distribution 
function /,-. The distribution function evolves according 
to the Boltzmann equation, 

fi(x + Vi,t + l)-fi(£,t) = Sli(3,t) , (1) 

where f2j (x, t) is a collision term, which is usually taken 
as a time relaxation to some equilibrium density, ff q . 
This is known as the the Bhatnagar-Gross-Krook (BGK) 
operator (l3j . 

n i {x,t) = --{f i {x,t)-ft\x,t)) , (2) 

T 

where r is the relaxation time and ff q (x, t) is the equilib- 
rium function. The equilibrium function is chosen in such 
a way, that (in the continuum limit) the model simulates 
the actual physics of the system. 

For our 3D model, we use a cubic regular grid of lat- 
tice constant Sx—\/2cSt, with c the light speed in vaccum 
(c~3 x 10 8 m/s); in other words, c=l/v / 2 in normalized 
lattice units (time unit = St, spatial unit = Sx). There 
are 13 velocity vectors (figure [TJ, 13 different vectors for 



FIG. 3: Index relationship between the velocity vectors and 
the electric and magnetic vectors. 

the electric field (figure [T]) and 7 different vectors for the 
magnetic field (figure ^ . The velocity vectors are de- 
noted by vf , where z=l, 2, 3, 4 indicates the direction and 
p=0, 1, 2 indicates the plane of location (figure [T]). They 
have magnitude y/2, in lattice units, and lie on the diag- 
onals of the planes. Their components are 

fl? = V2(cos((2i-l)7r/4),sin((2i-l)7r/4),Q) , (3a) 
vj = V2(cos((2i - l)7r/4), 0,sin((2i- 1)tt/4)) , (3b) 

v\ = y/2(0, cos((2i - 1)71-/4), sin((2i - 1)tt/4)) . (3c) 

This makes 12 vectors. The missing one is the rest vector 
vq, with components (0,0,0). 

Associated to each velocity vector vf there are two 

electric vectors and two magnetic vectors (j=0, 1), 
which are used to compute the electromagnetic fields (fig- 
ure [3]). The electric vectors are perpendicular to iff and 
lie on the same plane p. The magnetic vectors are per- 
pendicular to vf, too, but they are also perpendicular to 
the plane p. In terms of the velocity vectors ©, they 
are: 

-p 1 -p -p 1 -/p 

e io ~ 2 V l( i + 2 ) mod ' &il ~ 2, V ^ mod 4 1 +1 ' 

" (4) 

with the rest vector e*o=(0,0,0), and 

%=ifx^ , (5) 

with the rest vector &o=(0,0, 0). With these definitions, 
there are 25 electric field vectors, but only 13 of them are 
different. Similary, there are 25 magnetic field vectors, 
but only 7 are different. 

The electromagnetic fields are computed from distri- 
bution functions that propagate from cell to cell with 
the velocity vectors m. There are four distribution func- 
tions associated with each velocity vector, denoted by 
G%j (i=0, 1 and r=0, 1), plus two functions associated 

(r) 

with the rest vector Vq, denoted by G . That makes 
4 x 12 + 2=50 distribution functions. 
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The macroscopic variables are computed as follows: 

4 2 1 



^ = EEE G 

i=l p=0 j=0 



<3 ij 



4 2 1 
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i=l p=0 j=0 



(6a) 



(6b) 



and the mean density current vector J' is defined by 

J' = aE' . (11) 

In order to have the mean density vector J' in terms of 
the subsidiary fields, we replace Eg.(|10p into Eq. (fTT|) to 
obtain 



J' 



(12) 



i=l p=0 j=o 



(6d) 



Next, we adopt BGK collision terms fl?- and fig of 
the form [H 



ng r) = -;(^ (r) (^*)-Gg r)Bq (^*)) . (13a) 



« = -i( G W(f,t)-4^(f,t)) 



(13b) 



fir 



J = oE 



(6e) 



where D, B and J are subsidiary fields that represent the 
displacement field, the electric field and the total current 
density before external forcements, respectively (the ac- 
tual mean fields, including external forcements, are de- 
scribed below). In addition, B is the induction field, H 
is the magnetic field, p c is the total charge density and 
e r , fj, r and a are the dielectric constant, the permeability 
constant and the conductivity of the medium, repectively. 

For the evolution of these distribution functions we 
follow the proponsal of Zhaoli Guo, Chuguang Zheng and 
Baochang Shi [3| that includes external forcements in 
the lattice Boltzmann equations, as follows: 



nf\x,t)+T^\ 



(J) 



where the relaxation time is set to t=|. 

The equilibrium functions for the electromagnetic 
fields are 



(6f) G^ )eq (f,i) 



G^(x,t)= Pc , 



(14a) 
(14b) 



This completes the definition of the lattice Boltzmann 
model. 

The proof that this latticc-Boltzmann model, via a 
Chapman-Enskog expansion, recovers the Maxwell equa- 
tions is shown in Appendix [X] This model reproduces 
the following equations: 



dp c 
dt 



V • J' = 



dB 
"dt 



(15a) 



(15b) 



G^ (x,t+l)- G rj (x, t) = n%> (x, t) + T V 



-oW/ 



(8) 



where and Tq ' are forcement coefficients (r=0, 1), 
defined by Q 



(r) 



T, 



The mean electric field, E', is defined by 

4e r 



(9a) 



(9b) 



(10) 



V x H = p Q J' 



1 dD' 
~&~dt 



(15c) 



It is known [l5| that we can derive from these equations 
the following expresions for the divergence of the electro- 
magnetic fields: 



V • D' 



Pc 

eo 



V • B = 



(16a) 



(16b) 



that is, the other two Maxwell equations hold if they are 
satisfied as initial conditions. 
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FIG. 4: Electric field at t=0. The shadow zone is the dielec- 
tric medium which it has a dielectric constant e=2.5 and the 
normal zone correspond to vacuum zone (e=1.0). 
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FIG. 5: Electric field at £=140 clicks. The shadow zone is 
the dielectric medium which it has a dielectric constant e=2.5 
and the normal zone correspond to vacuum zone (e=1.0). 



III. APPLICATIONS 

A. Dielectric Interface 

First, we simulate the propagation of a electromag- 
netic gaussian pulse crossing an dielectric interface. For 
this purpose, we took an unidimensional array of L cells 
with periodic boundary conditions in the z coordinate 
and with each cell being its own neighbor in both direc- 
tions of the x and y coordinates. One half of the sim- 
ulation space, z<L/2, is vaccum (e=1.0) and the other 
half, z>L/2, represents a dielectric medium with dielec- 
tric constant e=2.5. For avoiding abrupts changes of 
the dielectric constant between two neigboring cells we 
choose the following distribution of the permitivity: 

e = 1.75 + 0.75 tanh(x - L/2) . (17) 

The functional form of the incident gaussian electromag- 
netic pulse is given by the equations 

E = (£ cxp(-aO-z ) 2 ),0,0) , (18a) 

B= (O,B o cxp(-a(0-z o ) 2 ),O) , (18b) 

centered at zq. The constant a fixes the pulse width, 
Eq is the pulse amplitude and the constant Bq is related 
with Eq according to the relation Eq=cBq, with c the 
vacuum light speed. 

For the simulation we choose £=200, c=l/v / 2, 
£7o=0.001, a=0.01 y zq=40 (in normalizated units). The 
initial condition is shown in figure^ and the electric field 
after 140 time steps is shown in figure O 

The theoretical predictions of the amplitudes of the 
transmited and reflected pulses can be calculated from 

m 

yt + 1 




(19b) 



For our case these predictions give a ratio between the 
transmited amplitude and the incident amplitude of 
0.7751 and a ratio between the reflected and incident 
amplitudes of 0.2249. The values measured in the sim- 
ulation are 0.7750 and 0.2249, respectively. This results 
correspond to errors smaller than 1% with the theoreti- 
cal results. The results of this simulation show us that 
the model is able to reproduce the propagation of elec- 
tromagnetic fields when dieletric media are present. This 
simulation takes less than 1 minute in a single Pentium 
IV de 3.0 GHz. 



B. Skin Effect 

For reproducing the skin effect, we simulate an elec- 
tromagnetic plane wave that arrives perpendicular to a 
conductor material. For this purpose we construct a uni- 
dimensional space of L cells as before, with zero conduc- 
tivity for z<L/4 and cto conductivity for z>L/4 is do- In 
order to avoid abrupt changes between neigbohring cells 
we choose the following distribution for the conductivity: 

cr = cr (l + tanh(ai - L/4)) . (20) 

In order to generate an incoming plane wave we set at 
z=0 an electric oscillator by imposing there an electro- 
magnetic field given by 

E= (E sm(u;t),0,0) , (21a) 

B = (0,£ sia(cj*),0) , (21b) 

where Eq is the wave amplitude, cBq—Eq and uj is the an- 
gular frequency of the wave. Some simulation snapshots 
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FIG. 6: Electric field E as a function of the z component with 
a conductor position z>250 at t=635 time steps. The solid 
fine is the electric field and the dashed line is the theoretical 
exponencial decreasement. 
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FIG. 7: Electric field E as a function of the z component with 
a conductor position 2>250 at t=647 time steps. The solid 
line is the electric field and the dashed line is the theoretical 
exponencial decreasement. 

are shown in Figures El [7] and [8] The simulation values 
are (in normalized units): L=1000, £Jo=0.001, w=7r/100, 
c=l/V2 and cr =10 6 . 

The theoretical expression for the amplitude of the os- 
cillating electric field inside the conductor is 

E = E exp(-x/A) , (22) 

where Eq is the amplitude of the electric field at the 
outside and A is known as the skin thickness. For good 
conductors this thickness is given by [l5| 

A= J—^— • (23) 

V afJ,Q/J, r UJ 

Figures El [7] and [5] also shows Eq. (|22p as a dashed line. 
One can observe that the amplitude of the electric field 
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FIG. 8: Electric field E as a function of the z component with 
a conductor position z>250 at t=670 time steps. The solid 
fine is the electric field and the dashed line is the theoretical 
exponencial decreasement. 

oscillation follows in excellent agreement the theoretical 
prediction. This simulation takes less than 1 minute in a 
single Pentium IV de 3.0 GHz. 

C. Electric Dipole 

To simulate an electric dipole we construct an array 
of 100 x 100 x 100 cells with free boundary conditions 
(each limit cell takes itself as his own missing neighbohr). 
In the center of this array we insert a small oscillating 
current density in the z direction, 

J z = JoSm(^jt\ , (24) 

where Jg is the amplitude of the current density and 
T is the oscillation. For avoiding any abrupt change of 
the physical quantities between two neigboring cells we 
choose actually a gaussian functional form for the ampli- 
tude of the current density J' , 

4 = J cxp(-0.75[(a; - 50) 2 + (y - 50) 2 + (z - 50) 2 ]) . 

(25) 

The periode was set to T=25.0 time steps and the am- 
plitude to J o =0.0001 (in automaton units). 

The simulation results are shown in Figures [H] and 1101 
Figure [9] draws the equipotential lines after 56 time steps. 
Figure [10] shows the amplitude of the radiated magnetic 
field as a function of the x component and the theoretical 
values given by [15j 

Here, k is the magnitude of the wave vector, k=^, with 
c the light speed in vaccum and u> the oscillation angu- 
lar frequency. The magnitude P of the electric dipolar 
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Simulational amplitude 
values B exp 


Theoretical amplitude 
values B theo 


Error 

(%) 


23.33 


23.01 


1.4 


9.53 


9.37 


1.7 


6.12 


6.07 


0.8 


4.54 


4.51 


0.7 



TABLE I: Simulational B exp and theoretical B t h eo magnetic 
field amplitudes for an oscillating electric dipole (all values in 
automaton units). 
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FIG. 9: Equipotential lines produced by an electric dipole in 
z after two whole oscillations (approximately). 



100x50x50 cells with free boundary conditions and we 
insert two parallel metallic layers of conductivity <r, one 
wider than the other. The dimensions, depicted in figure 
ITT1 are (in normalizated units): t=5, d~20, W=50 and 
e r =l. The signal enters into the waveguide by forcing 
35, ■ 1 




60 65 70 75 
Z (cells) 



85 90 



FIG. 10: Amplitude of the oscillating magnetic field produced 
by the electric dipole of figure [9] The solid line represents the 
theoretical value (Eq. ([26}) and the dashed line is the one 
obtained from the simulation. 




momentum is computed as 



P = 



(27) 



where £ is the effective volume of the dipole. Table U 
shows the simulation peaks and the theoretical values 
for the magnetic field amplitude. Again, the simula- 
tion matches the theoretical predictions with an accuracy 
of less than 2%. This illustrates the possibilities of the 
method to simulate antennas. This simulation takes 40 
minutes in a single Pentium IV de 3.0 GHz. 



D. Waveguide. Microstrip 

In this section we simulate a microstrip waveguide (see 
figure [TTj) . For this purpose, we chose a 3D array of 



FIG. 11: Microstrip waveguide. Here W is the width of the 
upper metallic plate and t its thickness, d is the dielectric 
thickness between the two parallel plates with permitivity e r . 



the cells at y—Q to have the electric and magnetic fields 
of a plane wave (Eqs. (|2"Tj) ) with i? =0.01, c=l/V2 and 
cj=7t/100 (see figure [T3|). The conductivity is smoothly 
changed along three cells, as in Sec. (see figure [T2"| . 

Figure [14] shows the z component of the electric field 
at y=25 after 467 time steps. The voltage and current 
along the waveguide at the same time are drawn in fig- 
ure [151 The electric voltage V(z,t) and current I(z,t) 
are in phase (plane wave). The waveguide impedance Z, 
computed as Z=V/I, gives us Z=70.73 fl The theo- 
retical value for the impedeance of an infinite microstrip 
waveguide is given by 
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X (cells) 

FIG. 12: Electric conductivity as a function of x and z at 
2/=25 




FIG. 13: Signal configuration that enter in the waveguide. 
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(28) 



where 



"e// = W + t- 



2tt 



In 



4611 1 



1 1 



t + 10 



(29) 



With this equation we obtain a theoretical value of 
Z=72.6f2. This gives us a 3% difference between the 
simulation and the theoretical predicition. This simula- 
tion takes 90 minutes in a single Pentium IV de 3.0 GHz. 



E. Resonant Cavity 

As a last test, we simulate a cubic resonant cavity and 
found its resonance frecuencies. The cavity we chose is 
an array of 50x50x50 cells, corresponding to a size of 
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FIG. 14: The z component of the electric field in the mi- 
crostrip at j/=25. 
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FIG. 15: Voltage (red line) and current (green line) in the 
microstrip waveguide at y=25 as a function of time. 



5x5x5 cm, with periodic boundary conditions, but im- 
posing a null electric field at the boundary (a perfect 
conductor). Inside the cavity we set an emmiter point 
as a single cell with an oscillating electric field in the x 
direction (the antenna) and a receptor point (the mea- 
surement instrument) as a sigle cell where we measure 
the electric field amplitude. The emiter point was set at 
(5, 5, 5) (in cell units) and the receptor point at (5, 5, 45). 
The oscillation frequency was scanned from 0.0182 to 
0.027 oscillations per time step by multipying each pre- 
vious value by a constant factor of 0.0072 (a total of 22 
different frequencies). Each frequency starts a single run, 
consisting of 5000 time steps before equilibrium and four 
whole oscillations to estimate the amplitude of the elec- 
tric field. Figure [16] shows amplitude of the electric 
field as a function of frequency. The resonant peaks are 
cleary identified. Table |TT] compares the values these res- 
onant frequencies with the theoretical predictions com- 
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Experimental values 


Theoretical values 


Error 


log(/exp) 


log (/the) 


(%) 


-1.725 


-1.724 


0.06 


-1.690 


-1.694 


0.24 


-1.660 


-1.669 


0.54 


-1.643 


-1.646 


0.18 


-1.618 


-1.625 


0.43 


-1.590 


-1.589 


0.57 



TABLE II: Simulational f exp and theoretical fthe resonant 
frequencies for a cavity of size 50x50x50 cells. The frequen- 
cies are in oscillations per time step. 



puted from (if 




pTT 



(30) 



All differences are smaller than 1%. The whole simu- 
lation (22 runs of more than 5000 time steps each on 
125000 cells) takes 18 days in a single Pentium IV of 3.0 
GHz (actually, we ran the simulation on a cluster of five 
nodes on two days). 



IV. CONCLUSION 

In this paper we introduce a three dimensional lattice- 
Boltzmann model to simulate the Maxwell equations in 
macroscopic media (dielectric, magnetic and contuctor 
materials). The model actually reproduces the Maxwell 
equations in the continuous limit, as prooved in Appendix 
lAl In order to test the model, we simulated some typical 
electromagnetic sets, namely: the amplitude of the trans- 
mitted and reflected pulses at a dielectric interface, the 
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exponential decay in the amplitude of a incident plane 
wave on a conductor (skin effect), the amplitude (as a 
function of distance) of the magnetic field radiated by 
an electric dipole, the characteristic impedance of a mi- 
crostrip waveguide and the resonant frequencies of a cu- 
bic resonant cavity. In all cases, the differences between 
the simulational results and the theoretical predictions 
lay between 1% and 3%. These five simulations proove 
that our LB model works very well in the more different 
situations. 

The model let us to select the electromagnetic con- 
stants (magnetic permeability, electric permitivity and 
conductivity) for each cell at pleasure. Nevertheless, for 
avoiding numerical instabilities on the interfaces between 
two different materials we need to smooth the transition 
across three cells, approximately. For this purpose we 
have used tanh functions, but other functional forms can 
be also employed. This is a standard procedure in many 
numerical models l2| . 

The main advantage of our LB model is that it is about 
ten times faster than FDTD based models. Indeed, a the 
simulation with FTDT of the propagation of an electo- 
magnetic pulse produced by a discharge antenna on a 
grid of tantos x tantas x tantas cells takes around Xh 
in a Pentium IV d. In comparison, our simulation of a 
Nos de que on a nose geu grid took xh. For this reason, 
the LB model seems to be promising in a variety of ap- 
plications, including among others: the propagation of 
electromagnetic pulses produced by atmospheric rays or 
discharge antennas, the light diffraction in objects that 
have a complex geometry or the sputtering of electromag- 
netic signals across buidings. The velocity and precision 
of our lattice-Boltzmann model for electrodynamics are 
excellent for the design and optimization of antennas, 
which requires to simulate a large number of configura- 
tions, and very specially in the design of pulse antennas, 
where single-frequency numerical methods fail. All these 
applications can be thema of future works. 

The developments in this paper also show how to intro- 
duce antisymmetrical tensors in conservation equations 
in Lattice-Boltzmann models. These equations are use- 
ful to simulate many interesting phenomena, like seismic 
or torsion waves. These are also topics of future devel- 
opment. 

Hereby we have introduced a Lattice-Boltzmann model 
for electrodynamics that actually reproduces the Maxwell 
equations in media, with a plenty of future aplications. 
Moreover, this model shows how to construct LB models 
that fulfills conservation laws with antisymmetric ten- 
sors. We hope that this valuable theoretical develop will 
push the evolution of LB models further away in the hori- 
zon of even more exciting applications. 



Boltzmann. We also thanks the scholarship program of 
the National University of Colombia for finnantial sup- 
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APPENDIX A: CHAPMAN-ENSKOG 
EXPANSION 

The Boltzmann equations for each fluid, Eq. ([7]) and 
([5]), determine the system evolution. This evolution rule 
gives in the continuum limit the macroscopic differential 
equations that the system satisfies. This is known as the 
Chapman-Enskog expansion. 

To develop the Chapman-Enskog expansion, we start 
by taking the Taylor expansion of these equations until 
second order in spatial and temporal variables, 



< :! ; ■ VG 



<p(r) 



1 d 2 G p 



dG 



<p(r) 
ij 



dt 
ld 2 G: 
2 



dt 



(Al) 



ij dt 2 



dt 2 



0G< r) 1 d 2 G 



(r) 



dt 



2 dt 2 



,p(r)eqN 
ij > 



(r)eqN 



(A2) 



where a, (3=x, y, z denote the x, y and z component. 

Next, we expand the distribution functions, the spa- 
tial and temporal derivates in a power series of a small 
parameter, A, 
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P(r) _ QP(r)(0) 



AG 



pM(i) 

ij 



2 n p(r)(2) 



AG 



— = \— A 2 — 

dt ~ dti + dt 2 



— - 9 

dx a dx a l 



(A3) 



(A4) 



(A5) 



It is assumed that only the Oth order terms in e of the 
distribution functions contribute to the macroscopic vari- 
ables. So, for ft > we have 
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The main current density J' is of the order A [TJ], so 



we can write J'=XJ'\. Because G 



■p(r)eq . 



is now function 



10 



of J', we need to develop a Chapman-Enskog expansion 
of the equilibrium function, too: 

QP(r)eq _ Qp(r)(0)eq \QP( r )W c( l _|_ )2QP(r)(2)cq (AJ) 



So, if we replace these results into Eqs. (|Alj) and (|A2|> . 
we obtain for the zeroth order in A 



r ,p(r)(0)cq _ n p(r)(0) 



(A8a) 



consider r=l/2. By summing up Eqs. (|A9ap , (|A9bp . 
(lAlOap and (|A10b|) over i, j and p, we get 



dp c 
dU 



= 



and 



dp c 



V • J'i = 



(A12) 



(A13) 



r> (r)eq _ ^(r)(0) 



(A8b) 



For the first order in A of the distribution function we 
have 



0f-V x G 



,p(r)(0) . dG 'ij 



p(r)(0) 



•pMC 1 ) _ QP(r)(l)eq 



(A9a) 



Now, we can add these two equations to obtain 
dp c 



dt 



V • J' = 



(A14) 



Multiplying the equations (|A9a|) . (|A9b|) . (|A10aP and 
(lAlObp by e?- and summing up over the index i, j and p, 
we obtain for r=0 

dUE) 1- /B\ 1 - f h \ 



AG, 



(r)(0) 



9ti 



= _I (G W(D„ G W(iM) , (A9b) 



and for the second order in A we obtain 
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(r)(0) 



dt 2 



(r)(2) 
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(r)(2)eq 




(AlOb) 



The first order and the second order terms for the 
equilibrium function of the electromagnetic fields are ob- 
tained by replacing the Eq. (fTtj|) into Eq. (fT4|) . From this, 
we can obtain, grouping in A powers, 



Gf m ° q (x,t) 



t E ■ e p - + - 

4 e « + 8 



B-% ,(Alla) 



and 



0(eE) 9 J'i 



dU 



4 9*i 



= . 



(A16) 



Summing up these two equations, and taking into ac- 
count the Eq. (fT0|) , we arrive to the first Maxwell equa- 
tion, 

d(eE') 1 - fB\ 1 - , k s 

J 9T i -2 VX (M) = - M °2 J ' (A17) 



Similary, multiplying the Eqs. (|A9ap and (|A10ajl by 
and summing up over i, j and p, we obtain for r=l 



and 



dB - 

— +Vx£=0 , 
9*i 



(A18) 



(A19) 



Summing up these two equations, we get the second 
Maxwell equation, 



Gf )(1)eq (x,*) = Afif . ^ - ^J'i • 4 ,(Allb) 



Gg r)(2)oq (f,*) = , 



G<p )c V,t) = Pc , 



(Allc) 



(Alld) 



G< r)(1)cq (f,*) = G^ )(2)oq (f,t) = . (Alle) 

Now, we are ready to determine the equations that 
the model satisfies in the continuum limit. First, let us 



dB 
~dt 



Vx£' = 



(A20) 



The others two Maxwell equations can be obtained from 
the Eqs. (|A17p and (|A20p how is shown in [T^]. If we 
applying the divergence to the Eqs. (|A19|) and (|A20[) we 

get 



9(V • E>) 1 



i)t 



MoV • J' , 



9(V ■ B) 

9* 



(A21) 



(A22) 
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and because the continuity equation Eq. (|A14[) . we arrive 
to 



<9(V ■ E') _ 1 dpc 
dt ~ 2 Mo ^t" 



Finally, we obtain 



0(V • E> - i M0 Pc) 



dt 



= o 



(A23) 



(A24) 



So, if the initial conditions of the electromagnetic fields 
satisfies the Maxwell equations 



V • B = 



(A25) 



^ 1 p r 

V • -E' = -MoPc = — 

l £q 



(A26) 



these equations will be reproduced correctly for all time. 



Summarying, the equations (|A17|) . (|A20j) . (IA25j) and 
(|A26[1 determine the electromagnetic fields evolution. 
These are the electrodynamics equations for macroscopic 
mediums. This complete the proof. 
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